Entregue os resultados na forma de um Jupyter Notebook, via Moodle, com todo o código, resultados (figuras e arquivos extras, se necessário), análise e conclusões.
Até o momento, nos preocupamos apenas com a magnitude da resposta em frequência de filtros seletores de sinais. Nesta atividade, iremos analisar o efeito da fase.
Considere um filtro passa-baixas analógico RC de 1a. ordem, com R=6,8 k$\Omega$ e C=22 nF.
from IPython.display import display, clear_output, Latex, Audio
import math
import cmath
import numpy as np
import scipy as scp
import scipy.io as sio
from scipy import signal
from lcapy import Circuit, jw
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.figsize':[15,9],
'font.size': 22,
'xaxis.labellocation': 'right',
'yaxis.labellocation': 'top'})
def rc_filter(R, C, Fs, filter_type):
# Returns equivalent IIR coefficients for an analog RC filter
#
# Usage: [B,A] = RC_FILTER(r, c, fs, type);
#
# R is the resistance value (in ohms)
# C is the capacitance value (in farrads)
# FS is the digital sample rate (in Hz)
# type is a character string defining filter type
# Choices are: 'high' or 'low'
#
# Highpass filters have the analog structure:
#
#
# | |
# Vi o--------| |----------+---------o Vo
# | | |
# C |
# ---
# | | R
# | |
# ---
# |
# |
# o---------------------+---------o
# GND
#
#
# Lowpass filters have the analog structure:
#
#
# |-----|
# Vi o--------| |------+---------o Vo
# |-----| |
# R |
# ----- C
# -----
# |
# |
# o---------------------+---------o
# GND
#
# This function uses a pre-calculated equation for both of these circuits
# that only requires the resistance and capacitance value to get a true
# digital filter equivalent to a basic analog filter.
#
# The math behind these equations is based off the basic bilinear transform
# technique that can be found in many DSP textbooks. The reference paper
# for this function was "Conversion of Analog to Digital Transfer
# Functions" by C. Sidney Burrus, page 6.
#
# This is also the equivalent of a 1st order butterworth with a cuttoff
# frequency of Fc = 1/(2*pi*R*C);
#
# Author: sparafucile17 07/01/02
#
# Verify that cutoff of the analog filter is below Nyquist
if( (1/(2*math.pi*R*C)) > (Fs/2) ):
raise ValueError('This analog filter cannot be realized with this sample rate');
# Default to allpass if invalid type is selected
b = [ 1, 0 ];
a = [ 1, 0 ];
# Constants
RC = R * C;
T = 1 / Fs;
# Analog Cutoff Fc
w = 1 / (RC);
# Prewarped coefficient for Bilinear transform
A = 1 / (math.tan((w*T) / 2));
####################################################
# The following equations were derived from
#
# s
# T(s) = -------
# s + 1
#
#
# using Bilinear transform of
#
# 1 ( 1 - z^-1 )
# s --> ----------- * ------------
# tan(w*T/2) ( 1 + z^-1 )
#
####################################################
if filter_type == 'high':
b[0] = (A) / (1 + A);
b[1] = -b[0];
a[1] = (1 - A) / (1 + A);
####################################################
# The following equations were derived from
#
# 1
# T(s) = -------
# s + 1
#
#
# using Bilinear transform of
#
# 1 ( 1 - z^-1 )
# s --> ----------- * ------------
# tan(w*T/2) ( 1 + z^-1 )
#
####################################################
if filter_type == 'low':
b[0] = (1) / (1 + A);
b[1] = b[0];
a[1] = (1 - A) / (1 + A);
return [b, a];
filtro = Circuit("""
P1 1 0; down=1.5, v=v_i(t)
R 1 2 6.8e3; right=1.5, scale=0.8, l={6,8\,\mathrm{k \Omega}}
C 2 0_2 22e-9; down, scale=0.8
W 0 0_2; right
W 2 3; right
W 0_2 0_3; right
P2 3 0_3; down, v^=v_o(t)
; draw_nodes=connections, label_ids=none, label_nodes=none""");
filtro.draw();
A. Qual a frequência de corte deste filtro?
H = filtro.P1.transfer('P2');
TH = H(jw);
display(Latex('$\\frac{V_{o}(s)}{V_{i}(s)} = H(s)$'));
display(Latex('$H(s) = $'+(H.canonical()).latex_math()));
display(Latex('$H(j\omega) = $'+(TH.timeconst()).latex_math()));
display(Latex('$\omega_{c} = \\frac{1250000}{187} = 6684,49197861 \, \mathrm{rad/s}$'));
display(Latex('$f_{c} = \\frac{\omega_{c}}{2\pi} \simeq 1063,86994045 \, \mathrm{Hz}$'));
B. Qual a magnitude da resposta em frequência exatamente na frequência de corte?
respcorte = TH.evaluate(1/(6.8e3*22e-9));
display(Latex('$|H(j\omega_{c})| = '+str(abs(respcorte)).replace('.',',')
+'\ \\left('+str(20*math.log10(abs(respcorte))).replace('.',',')+'\, \mathrm{dB} \\right)$'));
H.dB.bode_plot((1,10e4));
C. Qual a fase da resposta em frequência exatamente na frequência de corte?
fasecorte = cmath.phase(respcorte);
display(Latex('$\mathrm{arg} \, H(j\omega_{c}) = '
+str(round(fasecorte*180/math.pi))+'º$'));
H.bode_plot((1,10e4), plot_type='phase-degrees');
D. Monte esse filtro no VISIR e meça a magnitude na frequência de corte. Mostre as imagens do osciloscópio e explique como você chegou nos resultados.
entrada = 2.548;
saida = 1.808;
display(Latex('$|H(j\omega_{c})| = '+str(saida/entrada).replace('.',',')
+'\ \\left('+str(20*math.log10(saida/entrada)).replace('.',',')+'\, \mathrm{dB} \\right)$'));
A magnitude do filtro RC na frequência de corte é igual a 0,7095761381475667, resultado da equação $\frac{V_{o}(j \omega_{c})}{V_{i}(j \omega_{c})}$. Onde a tensão de pico de entrada, $V_{i}(j \omega_{c})$, é 2,548 V e a tensão de pico de saída, $V_{o}(j \omega_{c})$, é 1,808 V.
E. Ainda no VISIR, meça o atraso de fase na frequência de corte. Mostre as imagens do osciloscópio e explique como você chegou nos resultados.
$\mathrm{arg} \, H(j\omega_{c}) = -44,71º$
O atraso de fase do filtro RC na frequência de corte é igual a -44,71º. Esse valor pode ser obtido diretamente da função measure selecionado na opção phase no osciloscópio, no que o sinal de entrada foi escolhido na comparação com a saída. Outra maneira para verificar a fase é medir manualmente a diferença dos instantes de tempo entre o sinal de entrada e o de saída, e realizar a razão entre essa diferença e o período de qualquer sinal enquanto é comparado com a razão do atraso de fase e do ângulo de um ciclo completo.
F. Qual o atraso de fase para uma frequência 10 vezes menor que a frequência de corte?
$\mathrm{arg} \, H \left( \frac{j\omega_{c}}{10} \right) = -5,333º$
O atraso de fase do filtro RC na década antes da frequência de corte é igual a -5,333º. Similar ao item anterior, esse valor foi obtido diretamente da função measure selecionado na opção phase no osciloscópio, no que o sinal de entrada foi escolhido na comparação com a saída.
G. Qual o atraso de fase para uma frequência 10 vezes maior que a frequência de corte?
atraso = 22e-6;
atrasofase = -atraso*360/((6.8e3*22e-9*2*math.pi)/10);
display(Latex('$ \mathrm{arg} \, H (10 \cdot j \omega_{c}) = ' + str(atrasofase).replace('.',',')
+'º$'));
O atraso de fase obtido do filtro RC na década após a frequência de corte é igual a -81,84º. Similar aos itens anteriores, esse valor foi obtido diretamente da função measure selecionado na opção phase no osciloscópio, no que o sinal de entrada foi escolhido na comparação com a saída. No entanto devido a baixa magnitude do sinal de saída, o sinal é afetado pelo ruído presente no ambiente, o que modifica a sua medida de atraso da fase. Caso se faça a medida manual, em que foi encontrado uma diferença de 22 microssegundos, o valor de atraso da fase é de, aproximadamente, -84,26º.
Projete um filtro IIR equivalente ao analógico. Use a função rc_filter de https://www.dsprelated.com/showcode/199.php, considerando frequência de amostragem de 22050 Hz.
H. Quais coeficientes você obteve para esse filtro IIR?
I. A partir dos coeficientes, mostre os gráficos de magnitude e fase. Comente os resultados, comparando com os valores e medições do filtro analógico.
fc = 1250000/(2*math.pi*187);
fs = 22050;
b, a = rc_filter(6.8e3, 22e-9, fs, 'low');
w, resposta = signal.freqz(b, a, fs=fs, worN=2**13);
fase = [cmath.phase(valor)*180/math.pi for valor in resposta];
magnitudec = np.interp(fc, w, abs(resposta));
fasec = np.interp(fc, w, fase);
display(Latex('$b ='+str(b).replace(',',';').replace('.',',')+'$'));
display(Latex('$a ='+str(a).replace(',',';').replace('.',',')+'$'));
display(Latex('$\mathrm{Magnitude \ na \ frequência \ de \ corte} ='
+str(magnitudec).replace('.',',')
+' \ ('+str(20*math.log10(magnitudec)).replace('.',',')+' \, \mathrm{dB})'+'$'));
display(Latex('$\mathrm{Fase \ na \ frequência \ de \ corte} ='
+str(fasec).replace('.',',')+' \mathrm{º}$'));
plt.figure();
plt.plot(w, 20*np.log10(abs(resposta)));
plt.title(r'$Resposta \ em \ frequência \ do \ filtro \ passa \ baixas$');
plt.ylabel(r'$Magnitude \ [\mathrm{dB}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.xscale("log");
plt.grid();
plt.tight_layout();
plt.show();
plt.figure();
plt.plot(w, fase);
plt.title(r'$Resposta \ em \ frequência \ do \ filtro \ passa \ baixas$');
plt.ylabel(r'$Fase \ [\mathrm{º}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.xscale("log");
plt.grid();
plt.tight_layout();
plt.show();
Foram obtidos praticamente os mesmos valores de magnitude e de fase na frequência de corte mas a limitação da largura de banda superior do filtro digital afeta a sua resposta em frequência próxima a metade da frequência de amostragem quando comparada ao filtro analógico.
Considere a situação em que um sinal de radar é corrompido por ruído aditivo que ocupa uma banda de frequências distinta daquela do sinal. Para recuperar o sinal de interesse, é necessário aplicar um filtro. Entretanto, como há requisitos temporais a atender para a identificação do tempo de atraso a partir do sinal refletido, é importante que o filtro não distorça o sinal recebido. Isso será explorado nesta atividade.
Uma figura de mérito usada para verificar as distorções de fase de um sistema é o atraso de grupo, definido como
$\tau (\omega) = \mathrm{grd} [H(e^{j\omega})] = - \frac{\mathrm{d} [\arg H(e^{j\omega})]}{\mathrm{d}\omega}$
Verifique qual é o atraso de grupo do sistema atraso ideal ([1] p. 168, capítulo 5, Análise no domínio transformado de sistemas LIT, seção 5.1.1, Fase e atraso de grupo da resposta em frequência, figura em anexo).
Considerando um sistema de atraso ideal como em [1]:
\begin{equation} h[n] = \delta[n-n_{d}] \end{equation}Cuja resposta em frequência é:
\begin{equation} H(e^{j\omega}) = e^{-j\omega \cdot n_{d}} \end{equation}Com isso a sua resposta em frequência da sua fase é:
\begin{equation} \arg H(e^{j\omega}) = - \omega \cdot n_{d} \end{equation}E pela definição anterior de atraso de grupo:
\begin{align} \tau(\omega) &= - \frac{\mathrm{d} [\arg H(e^{j\omega})]}{\mathrm{d}\omega} \\ \tau(\omega) &= - \frac{\mathrm{d} [- \omega \cdot n_{d}]}{\mathrm{d}\omega} \nonumber \\ \tau(\omega) &= n_{d} \nonumber \end{align}A atividade seguinte foi adaptada de [2]. Carregue o arquivo de dados gdeldata.mat. No Python, use a função loadmat do módulo scipy.io (documentação).
Considere o sinal apresentado no vetor x1. Esse sinal será corrompido por ruído aditivo, mostrado, para fins didáticos, no vetor noise.
J. Observe a faixa de frequências que o sinal e o ruído ocupam. É possível projetar um filtro que recupere o sinal?
dados = sio.loadmat("gdeldata.mat");
x1 = dados['x1'][0];
ruido = dados['noise'][0];
x1ruido = x1 + ruido;
fs = 44100;
fftx1 = np.fft.fftshift(np.fft.fft(x1)/len(x1));
fftx1 = fftx1[int(len(x1)/2)::];
freqx1 = np.fft.fftshift(np.fft.fftfreq(len(x1), d=1/fs));
freqx1 = freqx1[int(len(x1)/2)::];
fftruido = np.fft.fftshift(np.fft.fft(ruido)/len(ruido));
fftruido = fftruido[int(len(ruido)/2)::];
freqruido = np.fft.fftshift(np.fft.fftfreq(len(ruido), d=1/fs));
freqruido = freqruido[int(len(ruido)/2)::];
fftx1ruido = np.fft.fftshift(np.fft.fft(x1ruido)/len(x1ruido));
fftx1ruido = fftx1ruido[int(len(x1ruido)/2)::];
freqx1ruido = np.fft.fftshift(np.fft.fftfreq(len(x1ruido), d=1/fs));
freqx1ruido = freqx1ruido[int(len(x1ruido)/2)::];
plt.figure();
plt.plot(freqx1, 20*np.log10(abs(fftx1)), 'b');
plt.title(r'$Resposta \ em \ frequência \ do \ sinal \ x1$');
plt.ylabel(r'$Magnitude \ [\mathrm{dB}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.grid();
plt.tight_layout();
plt.figure();
plt.plot(freqruido, 20*np.log10(abs(fftruido)), 'm');
plt.title(r'$Resposta \ em \ frequência \ do \ ruído$');
plt.ylabel(r'$Magnitude \ [\mathrm{dB}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.grid();
plt.tight_layout();
plt.figure();
plt.plot(freqx1ruido, 20*np.log10(abs(fftx1ruido)), 'k');
plt.title(r'$Resposta \ em \ frequência \ do \ sinal \ x1 \ e \ ruído \ aditivo$');
plt.ylabel(r'$Magnitude \ [\mathrm{dB}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.grid();
plt.tight_layout();
plt.show();
Sim, é possível projetar um filtro passa-baixas, já que a faixa de frequência ocupada pelo sinal e pelo ruído são diferentes.
Dois filtros foram projetados. O filtro $H_1(z)$, de resposta finita ao impulso, é representado por sua resposta ao impulso no vetor h. Já $H_2(z)$, um filtro IIR, é representado por dois vetores de coeficientes, a e b.
K. Confirme que os dois filtros apresentam resposta em magnitude que permite a recuperação do sinal de interesse. Sobreponha as respostas de cada filtro em um único gráfico.
dados = sio.loadmat("gdeldata.mat");
h = dados['h'][0];
b, a = dados['b'][0], dados['a'][0];
fs = 44100;
wfir, magfir = signal.freqz(h, fs=fs);
wiir, magiir = signal.freqz(b, a, fs=fs);
plt.figure();
plt.plot(wfir, 20 * np.log10(np.abs(magfir)), 'b', label="Filtro FIR");
plt.plot(wiir, 20 * np.log10(np.abs(magiir)), 'r', label="Filtro IIR");
plt.title(r'$Resposta \ em \ frequência \ dos \ filtros \ FIR \ e \ IIR$');
plt.ylabel(r'$Magnitude \ [\mathrm{dB}]$');
plt.xlabel(r'$Frequência \ [\mathrm{Hz}]$');
plt.legend(loc="upper right");
plt.grid();
Ambos são filtros passa-baixas capazes de recuperar o sinal de interesse.
Observe a resposta ao impulso do filtro FIR.
L. É simétrica? Que informação sobre a fase é possível estimar a partir dessa resposta ao impulso? Calcule e mostre o atraso de grupo desse filtro (no Python, com scipy.signal.group_delay).
dados = sio.loadmat("gdeldata.mat");
h = dados['h'][0];
fs = 44100;
plt.figure();
plt.stem(h);
plt.title(r'$Resposta \ ao \ impulso \ do \ filtro \ FIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
freqdelay, grpdelay = signal.group_delay((h, 1), fs=fs)
plt.figure();
plt.plot(freqdelay, grpdelay, 'b');
plt.title(r'$Atraso \ de \ grupo \ do \ filtro \ FIR$');
plt.ylabel(r'$Atraso$');
plt.xlabel(r'$Frequência \ por \ amostra \ [\mathrm{Hz/amostra}]$');
plt.ylim(0, 22);
plt.yticks([0,2,4,6,8,10,12,14,16,18,20,22]);
plt.grid();
A resposta ao impulso é simétrica, isso significa que a resposta em frequência da fase é linear. O atraso de grupo constante comprova isto.
M. Calcule e mostre o atraso de grupo do filtro IIR.
dados = sio.loadmat("gdeldata.mat");
b, a = dados['b'][0], dados['a'][0];
fs = 44100;
freqdelay, grpdelay = signal.group_delay((b, a), fs=fs);
plt.figure();
plt.plot(freqdelay, grpdelay, 'r');
plt.title(r'$Atraso \ de \ grupo \ do \ filtro \ IIR$');
plt.ylabel(r'$Atraso$');
plt.xlabel(r'$Frequência \ por \ amostra \ [\mathrm{Hz/amostra}]$');
plt.grid();
Filtre o sinal $x_1(n)$ usando cada um dos filtros.
N. Desenhe e analise os sinais de saída.
dados = sio.loadmat("gdeldata.mat");
x1 = dados['x1'][0];
ruido = dados['noise'][0];
x1ruido = x1 + ruido;
h = dados['h'][0];
b, a = dados['b'][0], dados['a'][0];
fir = signal.lfilter(h, 1, x1ruido);
iir = signal.lfilter(b, a, x1ruido);
plt.figure();
plt.subplot(211);
plt.plot(x1, 'b', fir, 'r');
plt.title(r'$Sinal \ original \ e \ filtrado \ pelo \ filtro \ FIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplot(212);
plt.plot(x1, 'b', iir, 'm');
plt.title(r'$Sinal \ original \ e \ filtrado \ pelo \ filtro \ IIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplots_adjust(hspace=0.5);
Os dois sinais filtrados apresentam atraso, porém o sinal filtrado por IIR tem um maior atraso e maiores distorções.
Filtre o sinal pulse com cada um dos filtros.
O. Desenhe e analise os sinais de saída.
dados = sio.loadmat("gdeldata.mat");
x1 = dados['x1'][0];
ruido = dados['noise'][0];
pulso = dados['pulse'][0];
h = dados['h'][0];
b, a = dados['b'][0], dados['a'][0];
fir = signal.lfilter(h, 1, pulso);
iir = signal.lfilter(b, a, pulso);
plt.figure();
plt.subplot(211);
plt.plot(pulso, 'b', fir, 'r');
plt.title(r'$Pulso \ original \ e \ filtrado \ pelo \ filtro \ FIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplot(212);
plt.plot(pulso, 'b', iir, 'm');
plt.title(r'$Pulso \ original \ e \ filtrado \ pelo \ filtro \ IIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplots_adjust(hspace=0.5);
Similar ao item anterior o filtro IIR apresenta não apenas atraso mas como distorções devido a sua não linearidade da fase, e no filtro FIR há apenas um atraso.
Filtre o sinal pnd_1 com cada um dos filtros.
P. Desenhe os sinais de saída e estime o tempo de atraso do pulso de cada um dos sinais resultantes. Explique qual o melhor filtro, FIR ou IIR, para essa aplicação.
dados = sio.loadmat("gdeldata.mat");
x1 = dados['x1'][0];
ruido = dados['noise'][0];
pnd = dados['pnd_1'][0];
h = dados['h'][0];
b, a = dados['b'][0], dados['a'][0];
fir = signal.lfilter(h, 1, pnd);
iir = signal.lfilter(b, a, pnd);
plt.figure();
plt.subplot(211);
plt.plot(fir, 'r');
plt.title(r'$Pulso \ filtrado \ pelo \ filtro \ FIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplot(212);
plt.plot(iir, 'm');
plt.title(r'$Pulso \ filtrado \ pelo \ filtro \ IIR$');
plt.ylabel(r'$Amplitude$');
plt.xlabel(r'$Amostras$');
plt.grid();
plt.subplots_adjust(hspace=0.5);
O pulso filtrado por FIR possui um atraso de 16 amostras, como visto nos itens anteriores. Já o pulso filtrado IIR possui um atraso estimado de 31 amostras obtidas a partir do gráfico acima. Visto que o filtro IIR possui grandes distorções e um maior atraso do sinal original, o filtro FIR é o mais adequado. Isso se deve à fase linear do filtro FIR.
[1] Oppenheim, A. V.; Schafer, R. W. Processamento em tempo discreto de sinais. 3. ed. São Paulo:Pearson, 2012.
[2] ECE 535 Digital Signal Processing: Matlab Assignment 3: Group Delay, Spring 2001. Disponível em: https://ece.gmu.edu/~kwage/courses/ece535/spr01/matlab/proj3.pdf
[3] ROHDE&SCHWARZ Pulsed RADAR signal generation and measurements : Educational Note 1MA234_0e. Disponível em: https://cdn.rohde-schwarz.com/pws/dl_downloads/dl_application/application_notes/1ma234/1MA234_0e_PulsedRadarEduNote.pdf
[4] WIKIPEDIA. Matched Filter: Examples. Disponível em: https://en.wikipedia.org/wiki/Matched_filter#Examples